Normalization and Batch Effect Correction after the new Peak calling
16 December 2021
1 Introduction
The new peak calling round applied on the previous notebook significantly increased the number of the features identified in our dataset. Therefore, we must need to repeat the standard downstream analysis, including data normalization, dimensionality reduction analysis and batch correction to account for this change in the number of features detected.
2 Pre-processing
2.1 Load packages
library(Seurat)
library(SeuratWrappers)
library(Signac)
library(harmony)
library(tidyverse)
library(plotly)
library(plyr)
set.seed(333)2.2 Parameters
cell_type = "GCBC"
color_palette <- c("#1CFFCE", "#90AD1C", "#C075A6", "#85660D",
"#5A5156", "#AA0DFE", "#F8A19F", "#F7E1A0",
"#1C8356", "#FEAF16", "#822E1C", "#C4451C",
"#1CBE4F", "#325A9B", "#F6222E", "#FE00FA",
"#FBE426", "#16FF32", "black", "#3283FE",
"#B00068", "#DEA0FD", "#B10DA1", "#E4E1E3",
"#90AD1C", "#FE00FA", "#85660D", "#3B00FB",
"#822E1C", "coral2", "#1CFFCE", "#1CBE4F",
"#3283FE", "#FBE426", "#F7E1A0", "#325A9B",
"#2ED9FF", "#B5EFB5", "#5A5156", "#DEA0FD",
"#FEAF16", "#683B79", "#B10DA1", "#1C7F93",
"#F8A19F", "dark orange", "#FEAF16", "#FBE426",
"Brown")
# Paths
path_to_obj <- paste0(
here::here("scATAC-seq/results/R_objects/level_4/"),
cell_type,
"/",
cell_type,
"_annotated_peak_calling_level_4.rds",
sep = ""
)
path_to_save <- paste0(
here::here("scATAC-seq/results/R_objects/level_4/"),
cell_type,
"/",
cell_type,
"_integration_peak_calling_level_4.rds",
sep = ""
)2.3 Load data
seurat <- readRDS(path_to_obj)
seurat## An object of class Seurat
## 142257 features across 21447 samples within 1 assay
## Active assay: peaks_redefined (142257 features, 0 variable features)
## 1 dimensional reduction calculated: umap3 Visualize UMAP without batch effect correction
# Normalization, dimensionality reduction
seurat <- seurat %>%
RunTFIDF() %>%
FindTopFeatures(min.cutoff = 10) %>%
RunSVD() %>%
RunUMAP(reduction = "lsi", dims = 2:30)
DepthCor(seurat)
DimPlot(seurat,
cols = color_palette,
group.by = "level_5",
pt.size = 0.1)
# Visualize UMAP's confounders
confounders <- c("library_name", "sex", "age_group", "hospital", "assay")
umaps_before_integration <- purrr::map(confounders, function(x) {
p <- DimPlot(seurat, group.by = x, pt.size = 0.1)
p
})
names(umaps_before_integration) <- confounders
print("UMAP colored by GEM:")## [1] "UMAP colored by GEM:"umaps_before_integration$library_name + NoLegend()
print("UMAP colored by sex, age group, cell hashing status, sampling center and assay:")## [1] "UMAP colored by sex, age group, cell hashing status, sampling center and assay:"umaps_before_integration[2:length(umaps_before_integration)]## $sex
##
## $age_group
##
## $hospital
##
## $assay
4 Run and visualize Harmony’s integration
seurat <- RunHarmony(
object = seurat,
dims = 2:30,
group.by.vars = 'gem_id',
reduction = 'lsi',
assay.use = 'peaks_redefined',
project.dim = FALSE,
max.iter.harmony = 20
)
# Non-linear dimension reduction and clustering
seurat <- RunUMAP(seurat, dims = 2:30, reduction = 'harmony')p1 <- DimPlot(seurat,
cols = color_palette,
group.by = "level_5",
pt.size = 0.1)p2 <- DimPlot(seurat,
cols = c("darkblue","darkblue","darkblue",
"darkblue","lightblue","lightblue",
"red","red","red","gray","gray"),
group.by = "level_5",
pt.size = 0.1)# Visualize UMAP's confounders
umaps_after_integration <- purrr::map(confounders, function(x) {
p <- DimPlot(seurat, group.by = x, pt.size = 0.1)
p
})
names(umaps_after_integration) <- confounders
print("UMAP colored by GEM:")## [1] "UMAP colored by GEM:"umaps_after_integration$library_name + NoLegend()
print("UMAP colored by sex, age group, cell hashing status, sampling center and assay:")## [1] "UMAP colored by sex, age group, cell hashing status, sampling center and assay:"umaps_after_integration[2:length(umaps_before_integration)]## $sex
##
## $age_group